home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_100 / 138_01 / rks.c < prev    next >
Text File  |  1985-03-09  |  5KB  |  235 lines
  1. /*
  2.     rk4n.c                created: 07-Nov-83
  3.     authors:            A. Skjellum, 
  4.                     M. Roberts
  5.  
  6.     updated:            14-Nov-83
  7.                     by MJR
  8.  
  9.     Copyright 1983, 1984 (c) California Institute of Technology.
  10.     All rights Reserved.  This program may be freely distributed
  11.     for all non-commercial purposes but may not be sold.
  12.  
  13.     Purpose:
  14.     integrate M first order differential equations
  15.  
  16.     y'[i] = f[i](t;y[j=1...M])
  17.  
  18.     M equations.
  19.  
  20.     algorithm:    see Burden, Faires, Reynolds, p. 239-240
  21.             also see p. 205
  22.  
  23.         1. interval t = [a,b]
  24.         2. choose N > 0 as as partition of interval (N steps)
  25.         3. define step size h = (b-a)/N.
  26.         4. Initial conditions: (denote w[i,j] as approxes to y's)
  27.  
  28.             w[i,0] = alpha[i]
  29.             means: ith w at time zero is set to initial value
  30.                    alpha[i].
  31.  
  32.         5. Computing the w[i,j+1] from w[i,j] is done as follows:
  33.  
  34.             loop over i = 1 to M
  35.  
  36.                 compute k1[i] = h*f[i](t,w[1,j],...,w[M,j])
  37.  
  38.             end of loop
  39.  
  40.             loop over i = 1 to M
  41.  
  42.                 compute k2[i] = h*f[i](t+h/2,w[1,j]+.5*k1[1],...,
  43.                         w[M,j] + .5*k1[M])
  44.  
  45.             end of loop
  46.  
  47.             loop over i = 1 to M
  48.  
  49.                 compute k3[i] =h*f[i](t+h/2,w[1,j]+.5*k2[1],...,
  50.                        w[M,j] + .5*k2[M])
  51.  
  52.             end of loop
  53.  
  54.             loop over i = 1 to M
  55.  
  56.                 compute k4[i] =h*f[i](t+h,w[1,j]+k3[1],...,
  57.                        w[M,j] + k3[M])
  58.  
  59.             end of loop
  60.  
  61.             loop over i = 1 to M
  62.  
  63.                 w[i,j+1] = w[i,j] +
  64.                        {k1[i] + 2*k2[i] + 2*k3[i] + k4[i]}/6
  65.                                                                               
  66.                         end of loop                       
  67.  
  68. */
  69.  
  70. double    (*rk_function)();
  71. double    (*rk_source)();
  72. double    (*rk_store)();
  73. double    (*rk_kuttas)[4];
  74. double    (*rk_comp[4])();        /* tells us how to form k's */
  75.  
  76. /* functions called indirectly by k_calc()  */
  77.  
  78. double rk_1(n,j,i)    /* to provide compatibility with calling */
  79. int n;
  80. int j;
  81. int i;
  82. {
  83.     return ((*rk_source)(j,n));
  84. }
  85.  
  86. double rk_23(n,k,j,i)  /*  N is in { 0...M } = argument number.
  87.                Since we have one argument to FN per equation,
  88.                N will indicate which we are currently being
  89.                 asked to provide.  Same goes for other RK_x
  90.                functions below.  */
  91. int n;
  92. int k;
  93. int j;
  94. int i;
  95. {
  96.     return((*rk_source)(j,n) + .5*rk_kuttas[n][k]);
  97. }
  98.  
  99. double rk_2(n,j,i)
  100. int n;
  101. int j;
  102. int i;
  103. {
  104.     return(rk_23(n,0,j,i));
  105. }
  106.  
  107. double rk_3(n,j,i)
  108. int n;
  109. int j;        
  110. int i;
  111. {
  112.     return(rk_23(n,1,j,i));
  113. }
  114.  
  115. double rk_4(n,j,i)
  116. int n;
  117. int j;
  118. int i;
  119. {
  120.     return((*rk_source)(j,n) + rk_kuttas[n][2]);
  121. }
  122.  
  123.  
  124. /*  Here's the integrator!!! */
  125.  
  126.  
  127. double rk4n(function,wsource,wstore,m,a,b,n,alpha,t,kuttas)
  128.  
  129. double (*function)();    /* pointer to function which returns deriv info */
  130. double (*wsource)();    /* source of w[i,j] values */
  131. double (*wstore)();    /* function which stores w[i,j] values for us */
  132. int m;            /* number of equations */
  133. double a;        /* start of interval */
  134. double b;        /* end   of interval */
  135. int n;            /* number of points */
  136. double alpha[];        /* array of initial values */
  137. double t[];        /* array where we store times */
  138. double kuttas[][4];    /* n x 4 kuttas (k1,k2,k3,k4 i=1,...m) */
  139. {
  140.  
  141.     register int i; /* looping variable */
  142.     register int j;    /* looping variables */
  143.     double time;
  144.     double h = (b-a)/(double)n;    /* step size */
  145.  
  146.     double rk_k2(),rk_k3(),rk_k4();    /* include for emphasis */
  147.  
  148.     rk_function = function;
  149.  
  150.     rk_kuttas   = kuttas;
  151.     rk_source   = wsource;
  152.     rk_store    = wstore;
  153.     
  154.     rk_comp[0]  = rk_1;
  155.     rk_comp[1]  = rk_2;
  156.     rk_comp[2]  = rk_3;
  157.     rk_comp[3]  = rk_4;
  158.     
  159.     /*  First assign initial values  */
  160.  
  161.     for (i = 0;i < m;i++)
  162.         (*rk_store)(0,i,alpha[i]);
  163.  
  164.     /*  Now loop through the necessary loop, calculating
  165.         each K value for each equation in each time step. */
  166.  
  167.  
  168.     for(j = 0;j < n;j++)    /* n time steps */
  169.     {
  170.         time = a + h*(double)j;  /* compute time */
  171.         t[j] = time;         /* also save it */
  172.  
  173.         rk4_1n(m,j,time,h);
  174.     }
  175. }
  176.  
  177. /* rk4_1n(): compute one solution step for n equations */
  178.  
  179. rk4_1n(m,j,time,h)
  180. int m;
  181. int j;                /* time step we're working on */
  182. double time;
  183. double h;
  184. {
  185.     register int k;    /* k calculation loop */
  186.     register int i;    /* m equations loop */
  187.     double value;    /* temporary */
  188.  
  189.     for(k=0;k < 4;k++)    /* k compute loop */
  190.         k_calc(m,k,j,time,h);
  191.  
  192.     for(i=0;i<m;i++)    /* compute new w[i,j]'s j fixed here */
  193.     {
  194.         value = (*rk_source)(j,i) + .166666667*(rk_kuttas[i][0] +
  195.             2.0*(rk_kuttas[i][1] + rk_kuttas[i][2])
  196.                 + rk_kuttas[i][3]);
  197.                 /* value for w[i,j] */
  198.         (*rk_store)(j+1,i,value); /* save this hard got number */
  199.  
  200.     } 
  201. }
  202.  
  203. /* k_calc(): compute rk coefficients for fixed j */
  204.  
  205. k_calc(m,k,j,time,h)
  206. int m;
  207. int k;
  208. int j;
  209. double time;
  210. double h;
  211. {
  212.     register int i;
  213.     double tval;
  214.  
  215.     switch(k)
  216.     {
  217.         case 0:
  218.             tval = time;
  219.             break;
  220.         case 1:
  221.         case 2:
  222.             tval = time + .5*h;
  223.             break;
  224.         case 3:
  225.             tval = time + h;
  226.             break;
  227.     }
  228.  
  229.     for(i=0;i<m;i++)        /*  used to be 1;<=m  */
  230.     {
  231.         rk_kuttas[i][k] = h*(*rk_function)(j,i,tval,rk_comp[k]);
  232.     }
  233. }
  234.  
  235.